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ABSTRACT 

As a test-bed for future investigations of directly imaged terrestrial exoplanets, we 
present the recovery of the surface components of the Earth from multi-band diurnal 
light curves obtained with the EPOXI spacecraft. We find that the presence and lon- 
gitudinal distribution of ocean, soil and vegetation are reasonably well reproduced by 
fitting the observed color variations with a simplified model composed of a priori known 
albedo spectra of ocean, soil, vegetation, snow and clouds. The effect of atmosphere, in- 
cluding clouds, on light scattered from surface components is modeled using a radiative 
transfer code. The required noise levels for future observations of exoplanets are also 
determined. Our model-dependent approach allows us to infer the presence of major 
elements of the planet (in the case of the Earth, clouds and ocean) with observations 
having S/N > 10 in most cases and with high confidence if S/N > 20. In addition, S/N 
> 100 enables us to detect the presence of components other than ocean and clouds in 
a fairly model-independent way. Degradation of our inversion procedure produced by 
cloud cover is also quantified. While cloud cover significantly dilutes the magnitude of 
color variations compared to the cloudless case, the pattern of color changes remains. 
Therefore, the possibility of investigating surface features through light curve fitting 
remains even for exoplanets with cloud cover similar to the Earth's. 

Subject headings: Earth - scattering - techniques: photometric 
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Introduction 



Exoplanet research is now a rapidly maturing branch of astronomy. In particular, a number 
of recently discovered planets have masses of order the Earth's and offer the possibility of study- 
ing other worlds much like our own. The initial analysis of the Kepler mi ssion data indicates 
that the small transiting planets are even more abundant than large ones (iBorucki et al.l I2010I . 
20 111 ) w hich is consistent wi th the statistical trend for planets discovered by the radial-velocity 
method (jHoward et al.l 120101 ). Moreover, candidates of nearly Ea rth-sized exoplanets within the 



habit able zones of their host stars have been already claimed (e.g. IVogt et al.ll2010l : iBorucki et al 
20 111 ). Nevertheless a mere detection of such a planet would be far from convincing evidence for 
habitability. 

Far richer information concerning the habitability would come from the atmospheric and sur- 
face properties of such planets, which will be probed by photometric and spectroscopic observation 
of planets in optical and near-IR bands. The light in the bands is dominated by starlight scattered 
by the planets, and the wavelength-dependent reflectivity and absorption features could provide 
direct indications of the planetary surface and atmospheric p roperties. This is why a number 
of on- going and future efforts including SE EDS ( Tamural 12009) . TM T (w ith a future instrumen t 
for direct im aging of Earth-like pla nets, see Matsuo &: Tamiiral 2010 ). 03 ( Savranskv et al. 2010), 
See-COAST (jSchneider et ahlbood ). SPICA fe.g. lEnva et ahlbood ) plan to develop instruments for 
direct observations of exoplanets. 

Until we have the first scattered light observation for an Earth-like exo planet, we may us e 
the Earth itself as a useful test-bed for future investigations as proposed by iFord et al.l (|200ll ). 
One approach is to extract spatially integrated scattered light from Earthshine. In addition to 
several clear absorption lines of biologically important molecules like H2O, O2, O3 and CH4(e.g. 
Pes Marais et al.l |2002| ) , a marginal signature of the red-edge due to ve getation (sharp increase 



of reflectivity at wavelength longer than 0.75//ni, e.g. ISeager et al.l 120051') has been reporte d (e.g 



Woolf et al.ll2002l : I Arnold et al.ll2002l : iMontafies-Rodrfguez et al 



amount consistent with simulations (e.g. iTinetti et al 



2006a 



200e 



Hamdani et al. 2006) in an 



b|; iMontafies-Rodrfguez et al.ll2006l ). 



Another approach is based on remote-sensing obser vation of the Earth , notably multi-band 
photometric data obtained with the EPOXI spacecraft. I Cowan et al.l (|2009l ) performed principle 
component analysis (PCA) in order to reconstruct surface features from the EPOXI data. Their 
PC A inversion identified two major eigen modes, the "red" and the "blue" components, which they 
roughly interpreted as land and ocean, respectively. They attempted to map the red component 
on the planetary surface using the diurnal variation, and found that it roughly reproduces the land 
distribution of the Earth along the longitudinal direction. We note here that their analysis neglects 
the effect of clouds on the diurnal light-curve since they assumed that the cloud cover fraction is 
constant at all time. 



Oaklev &: CashI (|2009l ) also attempted a longitudinal mapping of land and ocean components 
from a single-band, instead of multi-band, simulated photometry using the difference in reflectivity 
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between land and ocean. 



Fujii et al.l ()20ld ) (hereinafter Paper I) developed another reconstruction method of the plane- 
tary surface. In contrast to model-independent approaches such as PCA, Paper I approximates the 
planetary surface by a combination of four components: "ocean," "soil," "vegetation," and "snow" . 
The fractional areas of the four components were fit to reproduce the simulated light-curves of 
multi-band photometry of the cloudless Earth, and the method recovers the presence of ocean, soil 
and even vegetation reasonably well. 

The present paper extends the methodology and analysis of Paper I in various respects: 
1) we now include a "cloud" component both in computing simulated diurnal multi-band light- 
curves (forward procedure) and in reconstructing the fractional areas of the five components (in- 
verse procedure) , 2) the s imulated light-curves are based on the radiative transfer code rstar6b 



( Nakaiima Sz Tanakalll988l 1^l instead of employing Paper I's single scattering approximation, and 



3) the reconstruction procedure is applied to the EPOXI data instead of simulated light curves. 

The rest of the paper is organized as follows. Section [2] describes the technical details of the 
EPOXI observations on which the subsequent sections are based. Section [3] discusses our forward 
procedure for computing the simulated light curves using rstardb with input of land albedo data 
obtained with the MODIS (MODerate resolution Imaging Spectroradiometer) onboard the Earth 
Observing Satellites Terra and Aqua plus cloud data from MODIS onboard the Terra. Our inverse 
procedure, i.e. the reconstruction of the fractional areas of different surface components, is detailed 
in Section [H We also consider six different combinations of assumed soil, clouds and atmosphere 
properties, and discuss the resulting systematic uncertainty of the reconstruction. Section [5] uses 
this improved methodology to investigate the effect of clouds and the limits they impose on the 
interpretations of planetary light curves. Section [6] summarizes our main conclusions. Appendix lAl 
supplements the description of pre-process of the input data, and Appendix iBl derives the longitu- 
dinal map reconstructed from EPOXI data via our inverse procedure. 



2. EPOXI observation of the scattered light from the Earth 



One significant improvement over Paper I of the present work is the use of direct observations 
of the Earth's diurnal light curves, both as a validation of our forward procedure and as input for 
the inverse procedure. In particular we employ observations of the Earth that were obtained during 
NASA's EPOXI mission. The EPOXI mission reus ed the Deep Impact spa cecraft to carry out multi- 

EPOXI observations of 



band photometry of the Earth in optical/near-IR (jLivengood et al 



2011 



the Earth covered one whole day in each of March, May and June of 2008. The May data include a 
lunar transit of the Earth and thus introduce complexities that are beyond the scope of the present 
effort; we therefore do not use the May data. 



OpenCLASTR: |http://www.ccsr .u-tokyo.ac.jp/~clastr/ 1 
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Fig. 1. — Snapshots of the observational geometry every two hours during EPOXI observations of 
the Earth on 18-19 March 2008 (top) and on 4-5 June 2008 (bottom). The black solid line at hour 
7 indicates the orientation of the Earth's rotational axis on each day. 
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Fig. 2. — Solid: filters of the high-resolution visible-light instrument (HRIVIS) on the Deep Impact 
spacecraft, used for the EPOXI mission. The 'violet' filter (nominally 0.35//m) is short pass, where 
the detector and optical system components limit short wavelength sensitivity to > 0.365/im. The 
'IR' filter (nominally 0.95/im) is long pass, which is cut off at ~ 0.97^m by CCD response and 
optical system components. Dashed: CCD response function. 

Figure [J depicts the geometry of the Earth during the March and June 2008 observing runs. 
The spacecraft observed the Earth 76.7% illuminated at a phase angle of 57.7° in March, near 
the northern spring equinox, and 61.5% illuminated at a phase angle of 76.6° in June, near the 
northern summer solstice. The data effectively mimic observations of Earth-twins near quadrature 
in exoplanetary systems because the distance between the spacecraft and the Earth in March 
(~ 2.7 X lO'^km = 0.18AU) and in June (~ 5 x lO'^km = 0.34AU) is a few thousand times the 
Earth's radius and thus represents a good approximation to the geometry of an observation from 
infinite range. We use photometry obtained with the Deep Impact HRIVIS instrument at one-hour 
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Fig. 3. — Apparent albedo of scattered light from the Earth observed by EPOXI. Left: the pho- 
tometry of 7 EPOXI bands averaged over the March and June data. Error bars show the standard 
deviation of the hourly flux measurements in each wavelength band and are dominated by actual 
variability, not measurement errors. Right: light curves of 0.37 — 0.4 (purple line), 0.41 — 0.5 (blue 
line), 0.5 - 0.6 (light blue), 0.6 - 0.7 (green), 0.7 - 0.8 (red), 0.8 - 0.89 (brown), and 0.9 - 0.97/im 
(black) in the March data (solid) and June data (dashed). The observational uncertainty is less 
than 0.1 %, which is negligible in this figure. 



intervals in each of seven bands that span the visible range: 0.37 — 0.4^m, 0.41 — 0.5/im, 0.5 — 0.6/im, 
0.6 - 0.7/im, 0.7 - 0.8^m, 0.8 - 0.89/im, and 0.9 - 0.97/im. 

Images obtained with the HRI VIS camera are proc essed automatically for bias-subtraction, 
flat-fieldi ng, an absolute calib ration (jKlaasen et al.ll2008l ). Data acquired during the EPOXI mis- 



sion (e.g. iBallard et al.l l2010l ) indicate that the pixel-by-pixel flat-fielding correction has become 
inaccurate by about 1%, with random distribution across the detector array. The EPOXI Earth 
observations combine signal from thousands of pixels per frame, reducing the significance of random 
deviations to negligibility. 

Photometric signal is extracted from EPOXI Earth images by aperture photometry, summing 
over a circular extraction region centered on Earth's disc. Background signal is estimated from 
an annular region of similar surface area. The calibrated image units multiplied by the known 
solid angle per pixel yields collected intensity in units of W/m^/^um. The intensity is corrected for 
changes in the range from Earth to spacecraft, which varies slightly over the course of each 24- hour 
observing sequence and varies significantly between observations. Variations in Earth's heliocentric 
distance, while minor, are corrected similarly. 



Measurement uncertainties in the photometry are very small, of order 0.1% (jLivengood et al 



2011 



) , evaluated from the distribution of pixel values in the annular background used in the photom- 
etry. Photometry of the Moon, which did not alter its aspect during a set of EPOXI measurements 
in May 2008, yields an independent estimate of photometric uncertainty of about 0.05% after scal- 
ing for the relative diameter of Earth and Moon. Actual variability of the Earth's signal is a much 
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greater source of variability in the measurements, creating a range of measurements that is approx- 
imately ±4-8% of the mean signal measured in each filter. This range of variability is defined by 
the limits of the central 68th-percentile of the photometric sample population, yielding "one sigma" 
confidence limits on the expected value to be determined from a single independent measurement. 



The complete data set and data-reduction procedures are described by lLivengood et al.l (|201ll ). 
The HRIVIS filter functions are shown in Figure [2] along with the CCD response function. Optical 
components within the telescope and spectrograph further limit the extreme blue and red sensitivity 
of the shortest and longest wavelengths bands, respectively. 
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Fig. 4. — Color-color diagram of EPOXI data from March (solid, red) and June 2008 (dashed, 
green). Numbers indicate hours from the start of the observation, as indicated in Figured) Points 
represent the colors of some representative Earth-surface components. 

In Figure O we di splay the EPOX I data interpreted in terms of apparent albedo p* . We follow 
the definition of p* in lOiu et al.l (j2003l ) , as the ratio of total scattered intensity to the intensity of 
the incident solar flux scattered by a lossless Lambert sphere of the planet's radius, at the observed 
phase angle, 

1(A) 



p*{X) 



|F,i?2^(a)' 

(vr — a) cos a + sin a 



vr 



(1) 



(2) 



where /(A) is the scattered intensity as a function of the wavelength A, is the incident flux, Rp 
is the planetary radius, a is the phase angle, and ^(a) is the phase function for scattering by a 
Lambert sphere. Note that it is 3/2 times the geometric albedo. 
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The left panel of Figure [3] displays the photometry of 7 EPOXI bands averaged over both of 
the observed days, with the standard deviation of hourly measurements indicated as error bars. 
The variability in the measurements does not result from finite observational precision, but instead 
is the actual variability of the Earth's signal. As mentioned above, the observational uncertainty 
of the EPOXI measurements results in negligible noise level (< 0.1% with no exposure time longer 
than 100 ms) because the source (the Earth) is extremely bright. The right panel of Figure [3] shows 
the diurnal variation of each color band in March (solid) and in June (dashed); this figure illustrates 
both the common patterns of and difference between the March and June data. These data are 
qu alitatively consist ent with the expected photometric variability of the Earth, as first discussed 



bv lFord et all (120011 ). 



The shape variations between the lightcurves in different filters suggest looking more closely 
at the systematic diurnal color variability. We define color quantitatively as 

Cafe = -2.5 logio(^), (3) 

where is the reflectivity in band x, following the conventional astronomical practice. We select 
the 0.41 — 0.5 and 0.6 — 0.7//m bands to measure the variation in the blue region of the spectrum 
at short wavelengths and the 0.6 — 0.7 and 0.8 — 0.89/im bands to quantify the variation in the red 
region of the spectrum at long wavelengths. Figure H] displays the resulting color-color diagram of 
the diurnal light curves of the Earth for the two days. The colors of some representative components 
as seen through clear atmospheric layers are plotted as points in Figure H] for reference. 

Unsurprisingly, the average colors of the Earth are similar to the colors of clouds, but somewhat 
displaced due to the surface components. 



3. Forward procedure: simulating the diurnal scattered light-curve of the Earth 

Our goal is to develop a methodology to reconstruct the surface properties from the diurnal 
scattered light curves of exoplanets in habitable zones. However, before proceeding to this inversion, 
we choose to gain better insights into the problem and to validate our assumptions and models by 
carrying out forward calculations and comparing the results to the EPOXI data. In particular, we 
aim to determine the accuracy and limitation of the models of surfaces, clouds, and atmospheres 
as well as to validate our simulation code. 



Therefore, in this section we combine the models and empirical data of the surface and the cloud 
distribution so as to reproduce the diurnal light curves in multi-bands observed by EPOXI. In addi- 
tion to the forward calcul ation presented in Paper I, numerous o th er groups have simulated the scat- 



tered light of the Earth ( 


Ford et al. 


2001: 


Tinetti et al. 


2006a 


b; 


Montahes-Rodn'euez et al. 


2006; 


Palle et al. 


2008; 


Oaklev & Cash 


2OO9I: Arnold et al. 


2009; 


Douffhtv k Wolf 


2010; 


Robinson et al. 



201lh . 
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Basically the computation of the light-curves of the Earth requires the integration of the 
scattered light over the illuminated and visible surface. We have to divide the surface into two 
dimensional pixels and sum up the contribution from each pixel by approximating it locally as a 
flat plane. We use 2.0° x 2.0° latitude/longitude pixels as a compromise between computation time 
and accuracy. 

Unlike Paper I, which neglects multiple scattering and molecular absorptions, we incorporate 
both effects by using the hne-by-line radiative transfer code rstarog This code has been developed 
intensively as a tool for interpreting the remote-sensing data of the Earth. It incorporates the 
detailed modeling of optical properties of the atmosphere and clouds. 



3.1. Model assumptions in the forward procedure 

3.1.1. Atmosphere and clouds 

Scattering processes in the Earth's atmosphere are complicated. Since the albedos of clouds 
are generally much higher than those of the planetary surface components, the broad distribution of 
the cloud optical depths crucially affects the light curve of the Earth. Atmospheric pressure is the 
key parameter that determines the optical depth of Rayleigh scattering as well as the atmospheric 
features due to molecules, and the spectral features depend sensitively on the composition of the 
atmosphere. Both affect the broad-band photometry considered here. 

For simplicity we adopt a representative set of parameters called the "US standard" model 
atmosphere and assume that the model is valid everywhere. Basically the model provides the 
molecular composition and the temperature-pressure relation as a function of the altitude. However, 
the total scattered light is not sensitive to the details of those profiles. 

In addition to the homogeneous "US standard" model atmosphere, our model includes an 
empirical inhomogeneous cloud distribution. The spatial distribution of clouds is characterized by 
the cloud cover fraction / and the cloud optical thickness Tcm in each pixel. 

We use the cloud data set from Terra/MODIS Atmosphere Level 3 ProduclEI ( Dorothv et al 



20061 ). These data are daily global maps based on MODIS onboard the Earth Observing Satellite 
Terra and are given in 1.0° x 1.0° latitude/longitude pixels. Thus we merge them into 2.0° x 2.0° 
pixels to match the resolution that we use in the numerical code. The cloud data are not available for 
some locations, and we interpolate the data according to the procedure described in the Appendix 
|A]to fill these gaps. 



^http://www.ccsr. u-tokyo.ac.jp/ clastr/dl/rstar6b.html 

■^U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, D.C., 1976 
*http: / /ladsweb. nascom.nasa.gov/ 
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Although we take account of the spatial distribution of clouds, their optical properties are as- 
sumed to be everywhere the same. These optical properties depend on various parameters including 
the vertical profile, composition, phase (liquid or solid), radii and shapes of particles. For simplicity, 
we assume that all of the clouds are located at the altitude of 3.5-6.5 km and are composed of pure 
water; we neglect aerosols. The size distribution function of cloud particles is assumed to be 




2(}n{a)y 



(1.0 X lO-^cm < r < 3.0 X lO'^cm) 
(otherwise) 



(4) 



where C is the normalization factor. It also uses the size parameter r„ 
Clouds with these properties are hereinafter called "middle clouds". 



8.0x 10~^cm and a = 1.5. 



3.1.2. planetary surface components: land, ocean and snow 

The scattering model for the planetary surface is basically the same as adopted in Paper I. 
However, there are small differences in our treatment of land Bidirectional Reflection Distribution 
Functions (BRDF). These are described below. 

The scattering properties of land are described by the Rossi-Li model, a parametrized BRDF. 
While the model has three different terms (isotropic, geometric, and volume terms; see Paper I for 
details), we now approximate the planetary surface as Lambertian (isotropic scattering surface), 
and neglect the geometric and volume terms. At the geometry of the EPOXI observation, this 
approximation leads to < 5% differences in reflectivity compared to the light curves simulated with 
all three BRDF terms. 

We also follow Paper I in using the BRDF coefficients from "snow-free gap-filled MODIS BRDF 
Model Parameters"EI, which is a spatially and temporally averaged product derived from the 0.05° 
resolution BRDF/albedo data. Therefore we again merge the data into 2.0° x 2.0° pixels by simply 
averaging the BRDF coefficients within each pixel (from 40 x 40 data points). 



The BRDF for ocean pixels is described in Appendix B of Paper I (see also iNakajimal Il983l ) . 
The key parameter in this BRDF is the wind velocity 10 m above the surface, which we set to 4 m/s 
everywhere. Note that this model neglects the sub-surface scatterin g. In reality, the sub- surface 



scattering results in a slightly bluer ocean color than our model (e.g. iMcLinden et al.lll997l ) 



Finally we take account of the snow cover on the planetary surface, neglected in Paper I. The 
fraction of the snow-covered area is very small, but it may affect the total scattered light curve 
because of its high albedo. For this purpose, we use a data set of the monthly mean ice/snow cover 
fraction from International Satellite Cloud Climatology Project (isccpfl The data are provided 



'''http://modis.gsfc.nasa.gov/ 
®http://isccp. giss.nasa.gov/index. html 
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in 2.5° X 2.5° grids, and we interpolate them into 2.0° x 2.0° pixels according to the nearest-grid- 
point assignment. We regard those pixels with snow cover fraction greater than 0.5 as snow pixels, 
and then assign the albedo of "fine snow" from ASTER spectral library (shown in the left panel of 
Figure [7] below). 

3.2. Radiative transfer calculation with rstarGh 

The scattering models for land, ocean and snow with atmosphere and clouds described above, 
allows us to run the radiative transfer code rstardb, and then sum up all of the contributions from 
the relevant pixels according to the assumed geometry of the star-planet-observer system. 

Apart from the surface BRDF model, the radiance of a pixel is determined by five parameters; 
the solar zenith angle (^o) for the incident ray (assumed to come in parallel), the zenith angle of 
the observation (^i) for the scattered ray, the azimuthal angle between the incident ray and the 
observer ((/)), the optical thickness (rdd) and the cloud cover fraction (/) (Tabled]). 

Our surface classification proceeds as follows; first we identify pixels of ocean, and assign the 
corresponding BRDF model for ocean which includes strongly anisotropic specular reflection. Next 
we adopt an isotropic term only from the Rossi-Li BRDF for the remaining pixels (i.e., land) as 
described above. The coefficient for the isotropic term is obtained from MODIS. Finally pixels 
covered by snow are identified according to the ISCCP data, then the albedo for those pixels is 
replaced by that of "fine snow" (shown in the left panel of Figure [7] below). 

Then the total radiance of a pixel is given as the sum of the radiance of the surface through 
clear and cloudy portions of the atmosphere above that pixel: 

Radtotai(6'o,6'i,</',7"cid,a) = (1 - /) Rad(6'o, 6*1, 0, a) + / Rad(6'o, 6*1, (/>, Tcm, a). (5) 

Therefore one has to compute Jiad^Oo, 9i, (p, Tdd, o,) for each pixel by solving the radiative transfer 
equation in principle. In practice, however, a significant gain in the computation time is achieved by 
tabulating Rad(0O) ^ij ''"cid) o) and computing the radiance for each pixel via a linear interpolation 

Table 1: Grid system for radiance look-up table. 
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of the resulting 5-dimensional look-up table. The grids for the five parameters are shown in Table 
[TJ This procedure is repeated for each of the 7 EPOXI bands shown in Figure [2j 



3.3. Comparison of simulated light-curves with the EPOXI data 

In the preceding subsections, we have described our forward procedure to compute the diurnal 
light curves of the Earth. We now compare the resultant simulated curves with the data taken by 
EPOXI on March 18-19th and June 4-5th, 2008 so as to validate the procedure. 

The geometric configuration of the Sun, the Earth and the observer is fixed for each EPOXI 
observing run. We use the cloud data from the Earth Observing Satellite Terra for the corresponding 
days. The incident flux from the Sun to the Earth is approximated to be plane-parallel, and we 
adopt the distant-observer approximation at the location of EPOXI (~ 5 x 10''km away from the 
Earth). In addition we can safely neglect both the orbital motion of the Earth (~ l°/day in the 
orbital phase angle), and the motion of EPOXI spacecraft during each observing run. 

Since our land BRDF coefficients are defined according to the MODIS bands, they do not 
strictly match the EPOXI photometric bands. We assign the MODIS data to the EPOXI bands as 
shown in Table [2j Three EPOXI bands (0.3 — 0.4/im, 0.7 — 0.8//m, and 0.9 — 1.0/um) do not have 
direct counterparts in the MODIS bands, which may affect the results to some extent. 

Then we simulate the instantaneous light curve at one hour intervals. Since the EPOXI expo- 
sure times are < 6.15 x 10~^ sec, this is an excellent approximation. 

FigureElcompares the amplitudes of the total radiance between the EPOXI data (solid) and our 
simulations (dashed) averaged over 24 hours; strictly speaking, we average 24 simulated radiances 
sampled instantaneously at every hour in each banqj. The agreement is within 10% except for the 



^Figures E] and [S] are in the same format as Figures 3 and 4 of lRobinson et al.l (|201ll ') 



Table 2: EPOXI photometric bands and the corresponding MODIS bands. 



EPOXI band[^m] 


central wavelength [/_fm] 


MODIS band[/im] 


0.3 - 0.4 


0.35 


0.459-0.479 


0.4-0.5 


0.45 


0.459-0.479 


0.5-0.6 


0.55 


0.545-0.565 


0.6-0.7 


0.65 


0.620-0.670 


0.7-0.8 


0.75 
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0.8-0.9 


0.85 
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0.9 - 1.0 


0.95 


0.841-0.876 
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Fig. 5. — The 1-day average photometry of 7 EPOXI bands (sohd) and those from simulated hght 
curves (dashed). The March data is shown in red and June in blue. Note that the difference in 
radiance between the March and June data is due to difference in phase, i.e., the illuminated and 
visible area is bigger for March observation. 



0.3 — 0.4/im band. Figure [6] plots the normalized diurnal light curves, i.e., the radiances each hour 
in each band divided by its 1-day average (plotted in Figure [5]) ; solid and dashed curves indicate the 
EPOXI and simulated data, respectively. The observed and simulated diurnal patterns are in good 
agreement, mostly within 5% accuracy. Even in the three bands (0.3 — 0.4^m, 0.7 — 0.8/.im, and 
0.9 — 1.0/im) in which the EPOXI and MODIS band correspondence is not so good, the difference 
is within 10%. 

Thus we conclude that our forward procedure reproduces the observation adequately for our 
purposes. 

In addition to the forward procedure adopted here, we have investigated several alternatives. 
In particular we simulated light curves by treating all cloud pixels as (1) isotropic gray scatterers 
of albedo 0.4 which completely obscure the surface or (2) Tpih = 10 clouds. Instead of rstardb, we 



employed (1) an analytic 2-stream approximation (iMeador fc Weaver! Il980l ) or (2) an alternative 
numerical code libRadtrai^. Finally we calculated light curves (1) without account for the specular 
reflection of ocean water or (2) based on an alternative set of MODIS reflection spectra of surface 
components, called "white sky" albedo. None of these variations reproduced the EPOXI data more 
successfully than the adopted procedure and some gave considerably inferior performance. 



'http: / /www.libradtran.org/ 
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Fig. 6. — Comparison of EPOXI (solid) and simulated (dashed) light curves, each normalized by 
its time average. Left: March. Right: June. 
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4. Inverse procedure: recovering the planetary surface 



4.1. 



Method 



So far we have developed a forward procedure to compute the multi-band diurnal light curves 
of the Earth, and confirmed that our model reproduces the EPOXI data reasonably well. Next we 
consider an inverse procedure to recover the Earth's surface components from EPOXI multi-band 
photometric data. This exercise is intended to model the recovery of planetary surface information 
from future photometric data for terrestrial extrasolar planets. At this point, it is necessary to 
employ additional assumptions so that the inversion is possible. In particular, we assume that the 
planetary surface consists of a finite number of scattering components each with a priori known 
isotropic reflection spectra. Five components (ocean, snow, vegetation, soil and cloud) are included 
and used to fit the diurnal EPOXI photometric data by a linear combination of their isotropically 
averaged albedo spectra. 

Thus the apparent albedo p*j{ti) of the Earth at the j-th band at time ti is described as a linear 
combination of the fc-th components: 



where Djfc is the effective albedo of the k-th component at the j-th band. We define the effective 
albedo as the ratio of the total scattered intensity of a sphere wholly covered with a single component 
and the cloudless atmosphere to that of a lossless (i.e., non-absorbing) Lambert sphere at the 
same phase. The coefficient, Ak{ti), corresponds to the geometrically weighted area of the k-th 
components: 



In the above expression, 'dQ{6,(j)) is the solar zenith angle viewed from the planetary surface coor- 
dinate {9,(p), 'di{6,(p) is the zenith angle of the observer, mk{9,4)) denotes the fractional area of 
the k-th type at the coordinate {9,(j)), and the integral is performed over the illuminated (I) and 
visible (V) area of the Earth from the observer. 

We should note here that equation ([6|) treats the cloud cover simply as a single additional 
component to the other 4 surface components. In reality, the clouds are very complicated and 
their scattering spectra are dependent on their thickness, altitude, composition, and particle size 
distribution. In addition, the scattering properties of clouds and those of the surface components 
below interact in a complex way. Thus a fully realistic model would be underdetermined and 
therefore is impractical. Instead, we choose to characterize the clouds by a single representative 
refiection spectrum with a fixed set of properties and a single optical depth (our fiducial model 
is a middle cloud with r^d = 10; see Table U]). The validity of this major simplification will be 
examined in section [5^21 




(6) 




(7) 
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Our inverse procedure is carried out by adjusting the values of Ak{ti) in equation ([6]) to optimize 
the fit. In principle, the fit is performed by minimizing x^'- 



with aj{ti) being the appropriate observational errors. However, the errors in the EPOXI data are 
very small, much smaller than those likely to be obtained in any future observations of exoplanets. 
Thus we employ the following as a measure of the goodness of the fit: 

The magnitude of x^{ti) implies the required observational photometric accuracy to achieve such 
a recovery of the surface properties (see Figure [T3|) . 

We imposed the condition Ak{ti) > for all k and U in order to avoid unphysical models. We 
mi nimized in equat i on Q usin g the Bounded Variable Least-Squares Solver (BVLS) developed 
by Lawson Sz Hanson (jl974l . EoQsl n. The BVLS is designed to solve linear least square problems 
with bounded conditions on variables. 



4.2. Albedo models 

The design matrix which corresponds to the effective albedo in the j-th. band of the A;-th 
component is the essence of our model. The initial model employs the five components: "ocean," 
"snow," "vegetation," "soil," and "clouds." The corresponding spectra are computed with the 
radiative transfer code rstardb assuming the US standard atmosphere, which takes account of the 
effects of molecular absorption and the atmosphere's pressure-temperature profile. 



®The original code of the BVLS is available through NETLIB ( |http://www.netlib.org/lawson-hanson/index.html[ ). 



Table 3: Models adopted in the inversion procedure. 



model 


soil 


cloud altitude 


cloud optical thickness 


atmosphere 
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A 


middle 


10 


US standard 


II 
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middle 
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US standard 
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low 
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US standard 
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A 


high 


10 


US standard 
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A 


middle 
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US standard 


VI 


A 


middle 


10 


no O3 
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Fig. 7. — The effective albedo spectra for fine snow, vegetation, soil A ("class: inceptisol" ) , soil 
B ("class: entisol") and ocean. Left: Spectra from the ASTER library without atmosphere (solid 
lines). The reflectivity of oceans (including the scatter ing beneath the surfac e) is plotted as a 
dashed line for reference (courtesy of G. Tinetti; see also iMcLinden et alj (jl997l )). Right: Spectra 
with the US standard atmosphere (but neglecting the effect of cloud) at the phase of June 4th, 
2008 (phase angle ~ 76.6°) computed with rstarGh. Filled circles indicate the averages over the 
seven EPOXI photometric bands. 



Table 4: Assumed properties of three classes of clouds. 

low cloud middle cloud high cloud 
altitude [km] ^^3^ 3.5-6.5 6.5-20.5 
pressure [mbar] 680-1000 440-680 40-440 
component pure water pure water pure ice 

size parameter [cm] 8.4 x 10"'' 8.0 x 10"'' 2.0 x 10"^ 
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Fig. 8. — Effective albedo of different cloud types at the June phase (76.6°) calculated by rstardb. 
Red: low clouds of optical thickness 10, blue: high clouds of optical thickness 10, and black: middle 
cloud of optical thickness 1 , 5 and 10. The definitions for "low" , "middle" and "high" clouds are 
given in Table [H 

The left panel of Figure [7] plots the input spectra of the surface albedo (without atmosphere) for 
snow, vegetation and soil taken from the ASTER spectral library. We consider two representative 
models for soil (A and B) to investigate the effects of the variation in soil composition on the results. 
More specifically, we adopt "fine snow" for snow, "grass" for vegetation, "Dark yellowish brown 
micaceous loam (class: inceptisol)" for soil A, and "Brown to dark brown sand (class: entisol)" 
for soil B, where the nomenclatures quoted are those of the ASTER spectral library. The "soil" 
assumed in Paper I corresponds to soil B in this paper. The scattering model for ocean surface 
is described in Section 13.1.21 The right panel of Figure [7] is the effective albedo spectra including 
the effects of atmosphere given by rstardb at the phase of the June observing run. Filled circles 
indicate the averages over the seven EPOXI photometric bands for the incident solar spectra; these 
are the elements of the design matrix Dji^. 

The effective albedos for clouds are computed by considering cloud layers above a zero-albedo 
surface. This assumption is valid for large cloud optical depth Tdd, but underestimates the real 
albedo for smaller Tdd- Nevertheless this is a practical necessity to make the inversion feasible. 
Figure [8] shows the reflection spectra of different cloud optical thicknesses and altitudes. We define 
"low, " "middle," and "high" clouds as summarized in Tabled! the other parameters for the clouds 
are the same as adopted in the forward procedure. 

Figure [8] indicates that the cloud albedos are fairly insensitive to the wavelength. The cloud 
altitude primarily changes the strength of the absorption lines, and the cloud optical thickness 
changes the overall amplitude of the albedo, though not strictly in a linear fashion. Again, the 
filled circles indicate the elements of the design matrix Djk- 

Table [3] lists six models corresponding to different combinations of cloud and soil properties. 




First we consider our fiducial model (model I) in ^4.31 and then compare the results to the other 
models in ^4.51 



4.3. Decomposing the EPOXI data with model I 

We apply our fiducial model (model I in Table [3]) to the multi-band photometric data of the 
Earth observed by EPOXI, and infer the fractional areas of the five components. The results are 
plotted in Figure for the 1-day average data, and in Figure [TU] at 1-hour intervals over the diurnal 
cycle, respectively. In both figures, the left panel shows the result for the observation on March 
18-19th, 2008, while the right panel on June 4-5th, 2008. 

The filled circles in both figures indicate the resulting weighted area fraction Ak- The solid 
reference lines for "ocean," "soil," "vegetation," and "snow" are computed from the IGBP classi- 
fication map and ISCCP snow/ice cover data; we merge the original IGBP classification with 17 
classes into our 4 surface types by the same scheme as Table 2 of Paper I. Thus they correspond 
to the cloudless Earth. More relevant estimates are given by the dashed reference lines in which 
the surface pixels covered by the clouds (tcm > 10 for long dashed, and Tdd > 1 for short dashed) 
are removed from those components and regarded as cloud pixels (the case with r^id > 20 is also 
plotted in dotted lines in the "cloud" panels of Figure [T0|) . The differences in reference lines of the 
left and right panels arise from the different phase of the illuminated and visible part of the Earth 
(Fig[T]), different cloud cover, and seasonal snow cover on the two different observing runs. 
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Fig. 9. — Estimated fractional areas of the five surface components from the 1-day averaged data 
on March 18-19th, 2008 (left) and on June 4-5th, 2008 (right) for Model I. Filled circles and error 
bars are computed from 1000 Monte-Carlo realizations with S/N=20 (solid) and S/N=10 (dot- 
dashed). Dashed lines indicate reference values based on the combination of IGBP classification, 
ISCCP ice/snow cover data and MODIS cloud cover on each observing day. These are intended to 
represent the ground truth. The histogram shown by solid lines displays the cloud-free case. The 
long- and short-dashed histograms exclude surface pixels with T^id > 10 and Tcm > 1 clouds. 



This inverse procedure seems to work reasonably well. In particular, "ocean" and "clouds" 
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Fig. 10. — The same as Figure [9] but for the diurnal hght curve data. The fihed circles and the 
errors are computed from 1000 Monte-Carlo realizations with S/N=20 at every data point {i.e. 
each band and each time). 



are recovered as dominant compared to the other components, and the estimated fraction of each 
component is slightly but systematically smaller than the IGBP values (solid lines). 

Moreover, the diurnal variation pattern in Figure [10] qualitatively follows that of the actual 
surface components. For instance, the peaks of ocean at t=4 and t=19 correspond to the Pacific 
and Atlantic oceans, respectively. The peak of soil corresponds to the Sahara desert, and three 
peaks of vegetation to the Asian, African, and American continents. 

The variation pattern of clouds is complicated by changes in the distribution of cloud optical 
thickness with time. For example, the inversion reproduces the Tdd > 1 cloud curves best in 
amplitude, but matches Tdd > 20 cloud variation better in shape. 

We note in passing that the time variation pattern of the coefficient of the "blue" eigen spec- 
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trum extracted by PCA, interpreted bv lCowan et al.l (|2009l ) to roughly correspond to ocean, does 
not reproduce the global distribution satisfactorily. They hypothesized that this results from the 
time variation of clouds, which they assume to be time independent. Instead, the analysis pre- 
sented here separates the ocean and clouds a priori and takes account of the time variation of the 
cloud cover. This is probably why our inversion reproduces the variation pattern of ocean more 
successfully. This i mplies that the blue component extracted from the model- independent PCA by 
Cowan et al. (j2009l ) does not correspond to ocean alone, but to a combination of ocean and clouds. 



The time variati on of each component can be translated into its longitudinal distribution 
(jCowan (fe Agoll I2OO8I ) . By the proper combination of the orbital m otion of the planet s , it w ould 
be even possible to obtain a 2-dimensional map as first proposed bv lKawahara &: Fuiiil (|20ld ). In 
Appendix |B] we exhibit the longitudinal distribution of each component obtained from the diurnal 
variation patterns. 



4.4. Signal-to-Noise requirement 

In order to estimate the feasibility of the inversion procedure given the large observational 
errors anticipated in future observations of Earth-like exoplanets, we create mock data sets with 
signal-to-noise ratios of 20 and 10 using the EPOXI data. The variance of the estimated values is 
computed from 1000 Monte-Carlo realizations of the input spectra and is plotted as error bars in 
Figure [9l Similarly the plotted error bars in Figure [TO] assume S/N=20. 

The errors in the estimated cloud cover are fairly small, while those of ocean fraction are 
significantly larger. This is simply due to the difference in their reflectivity; the lower the reflectivity 
is, the less accurate the inversion becomes. The errors of vegetation are relatively small, probably 
due to the strong signature of the red edge. At least in the case of the Earth, the presence of clouds 
and ocean is successfully inferred from the time-averaged data with S/N ~ 10. However, their 
detectability is highly variable as the planet rotates and different portions of its surface are visible. 
So if the area of ocean on an extrasolar planet is smaller, its detection is therefore more difficult. 
For instance, a convincing detection of the presence of ocean with fractional area fraction of ~ 0.4 
(roughly corresponding to the Earth at t ~ 14 in the March and June observing run) would require 
S/N more than 20. 

We may want to estimate the minimum requirement for direct imaging instruments to attain 
this S/N. In an extremely idealized case where the observational uncertainty comes from the photon 
noise of planetary light only, the S/N is app roximately propo rtional to -^photon where photon 



counts A^'photon is scaled with equation (16) of lFujii et al.l (j20ld ). In this case, a 4 m telescope will 



be able to observe an Earth-twin at 10 pc away with S/N ~ 20 in 1 hour. The same S/N is also 
available with a smaller telescope if we carry out hourly exposures for several days and p hase them 



according to the planetary spin rotation period once it is determined (jPalle et al.ll2008l ) 



The photon counts will, however, substantially decrease through realistic future direct imaging 
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instruments, due to the throughput and the quantum efficiency of CCD. Moreover, the dominant 
noise hkely comes from speckles, exo-zodiacal light, and other instrumental noises such as readout 
noise and dark current. To be somewhat reahstic, we consider the effect of throughput, quantum 
efficiency, readout noise, dark current, and exo-zodiacal light. The assumed inst r umen tal design 
is listed in Table O which is basically the same as Table 1 of iKawahara fc Fuiiil (|20ld ) . On the 
assumption of these properties, S/N ~ 20 roughly corresponds to 4 m telescope observations of 
an Earth-twin at 10 pc away with 1-hour exposure with the data stacked for 6 days and phased 
according to the rotation period (evaluated at 0.8-0.9 /xm). The speckle noise depends on the details 
of the system which occults the light from the host star (coronagraph or external occulter). Since 
many such designs are under development, we do not try to account for noise due to residual light 
from the host star at this stage. Instead, the S/N requirement derived here will provide a goal of 
instrumental designs. 



Table 5: Assumed observational parameters 



Symbol 


Quantity 


Value 


Unit 


^' 


sharpness 


0.0433 




a 


pixel scale 


0.03125 


arcsec/pixel 


h 


end-to-end efficiency 


0.5 




V 


dark rate 


0.001 


counts/sec 


K 


read noise 


2 


"v/counts/read 


QE 


quantum efficiency 


0.91 






zodiacal light in magnitude 


23 


mag/square arcsec 


z 


zodiacal light 


1 





4.5. Inversion Sensitivity to the model assumptions 

As described earlier, our inversion may depend heavily on the adopted albedos for the selected 
surface components. It will not be easy to evaluate such systematic uncertainties associated with 
the model assumptions for future studies of exoplanets. However, we may gain some insights into 
the problem by comparing inversions of the EPOXI data using the different models listed in Table 

El 

Figure [11] is the same as the left panel of Figure [H but for models I to VI. Now the symbols 
indicate the estimates returned by the different models, and the vertical dashed lines indicate the 
range of model-to-model variations. 

Most of the model dependence can be qualitatively understood. Model II (cross), in which 
soil B is assumed instead of soil A, overestimates the area for soil because the effective albedo of 
soil B shows a lower-amplitude red portion of its spectrum (Fig[7|). Models III (asterisk) and IV 
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Fig. 11. — Decomposition of the averaged 7-band photometry of March with models I, II, III, IV, 
V and VI. The description of these models is given in Table [3l Histograms are the estimates of 
the ground truth identical to those in Figure [H Vertical dashed lines connect the maximum and 
minimum values returned by the various models. The statistical variations produced by signal-to- 
noise ratio 10 are shown as vertical solid lines. 



(open square) assume different altitudes for cloud layer. The difference in cloud altitude changes 
the absorption depths (Fig [8]) , and therefore the estimated areas of all the components are affected 
accordingly. If we adopt model V (plus) with Tdd = 5 clouds, which have a reflectivity smaller than 
those with r^d = 10, we end up with a larger estimate of the cloud area. Again, the estimation of 
all of the surface components is also affected by this change. Model VI (open circle) assumes no 
ozone in the atmosphere, and the lack of ozone's broad absorption at A ~ 0.55/im decreases the 
slope of all albedo models at short wavelengths. As a result, the area of ocean is estimated to be 
larger to compensate it. 

For reference we also plot the model I result with S/N=10 in Figure [TTJ The corresponding 
error bars for the fractional area are similar to the variation among different models considered here. 
This implies that the model-dependent systematic errors are not dominant for an observational S /N 
of order 10. 

This conclusion is also consistent with the result for the diurnal light curves shown in Figure [T2l 
The shaded region in each panel of this Figure indicates the range of the estimates among different 
models. Again the basic features are similar to those shown in Figure \W\ while the amplitudes of 
the estimated area vary fairly significantly from model to model, the overall variation patterns are 
relatively stable. Two of the minor surface components, "soil" and "snow", are difficult to identify 
robustly, but the major components, "cloud" and "ocean", are detectable in all the models. Finally 
the presence of "vegetation" may be detectable at some phases without being confused with the 
other components, probably due to its conspicuous "red edge" feature. 
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Fig. 12. — The range of decompositions of the June hght curves produced by models I, II, III, IV, 
V and VI. The shadowed region is bounded by the maximum and minimum values determined by 
this set of models. 



So far we have shown that our inversion procedure can identify ocean and clouds reasonably 
well but that the other components are not securely identified. Nevertheless, we may ask if an 
additional, unknown surface component is required to achieve a satisfactory fit. This is a retreat 
in a sense but may be a less model-dependent approach. 

For this purpose, we compute and compare the goodness-of-fit to the observed data with and 
without specific surface types. Figure [13] shows xi^i) per square-root of the degree-of-freedom for 
various models (without vegetation for simplicity). The reference model plotted in black solid 
line assumes that the planetary surface comprises of clouds, ocean and soil A. The other curves 
correspond to the four variants of the reference model: models with clouds, ocean and soil B instead 
of soil A (dot-dashed), with cloud and soil A without ocean (dotted), with ocean and soil A without 
clouds (long-dotted), and with clouds and ocean without soil (short-dotted). 
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Fig. 13. — The x of different models as a function of time. The curve labeled "reference" assumes 
ocean, soil A and middle clouds of optical thickness 10 but no vegetation component. The curve 
labeled "ref— clouds" means that the model omits clouds but otherwise is identical to the reference 
model, and the other labels are to be understood in the same manner. 

The resulting residuals of x/V d.o.f may be roughly interpreted as the inverse of the signal-to- 
noise ratio required to detect the discrepancy between the model and the observation. Since Figure 
[13] indicates that models with both ocean and clouds have values as low as 0.01, the presence of 
an additional component would be revealed by the data only with S/N > 100. On the other hand, 
models without either clouds or ocean show the difference in xl Vd.o.f values of 0.02-0.03 at certain 
phases of the light curves. It provides a rough idea of the signal-to-noise ratio (or the amount of 
the exposure time) required to identify such surface components in future exoplanet observations. 

Incidentally it is interesting to note that the time variation pattern of x/'V^d.o.f essentially 
follows that of the missing component; for instance, the model without ocean shows the largest 
discrepancy at i ~ 2 when almost all of the illuminated and visible area is covered by ocean (Fig{T]) . 
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5. Effects of clouds on photometric characterization of Earth-like Exoplanets 

Our primary goal in developing an improved treatment of forward and inverse procedures is 
to allow us to examine the effects of clouds. 

The presence of clouds affects the recoverable fractional areas of the different surface compo- 
nents in at least three ways: blocking the direct reflection spectra from the underlying surface areas, 
the nonlinear modification of the underlying surface reflection spectra by semi-transparent clouds, 
and the additional and strong time variability they introduce. We examine these three effects in 
turn below. 
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5.1. Cloud cover of the underlying planetary surface 



Approximately half of the surface of the Earth is covered by clouds with Tcm > 1. To the extent 
that the corresponding regions are completely unseen, the total surface area that one can directly 
observe is reduced by the same amount. Moreover the relative fractions of surface components 
contributing the light curves are m odified because ocea n and vegetation surfaces are more likely to 
be cloud-covered than soil (see also lCowan et alJl2009l ). 
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Fig. 14. — Color-color diagram. The trajectories of the simulated light curve of our "forward" 
model with realistic cloud cover, that of the "reduced cloud" model, and that of "no cloud model" 
are displayed together along with the representative colors of various components indicated by 
points. 

This effect of clouds may be evaluated more quantitatively with the forward calculation. To 
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be specific, we consider hypothetical Earth-twins with different cloud cover fractions; cloudless, 
reduced cloudiness, and realistic cloudiness cases. The realistic cloud case is based on the MODIS 
cloud data as described in Section [3l while the reduced cloudiness case has roughly 50 % of the 
cloud cover fraction. The latter is created by multiplying both the observed cloud cover fraction 
and cloud optical thickness at each pixel by a uniformly distributed random number between and 
1. The cloudless case neglects clouds entirely. 

The resultant color-color plot for 0.65/im and 0.85^m versus 0.45^m and 0.65/im is shown 
in Figure [TH The color trajectories over a diurnal period of the Earth are shown by the three 
curves. For reference, the colors of hypothetical planets with uniform surfaces of ocean (plus), 
soil A (open square), soil B (painted square), vegetation (asterisk), and snow (cross) with the US 
standard atmosphere are plotted as points. This figure indicates clearly why clouds dramatically 
limit our ability to recover the surface components from light curves. As the cloud cover fraction 
increases, the overall color approaches that of the clouds, and the amplitude of color variations is 
significantly reduced. The qualitative behavior is similar for color-color plots in other bands. 

It is, however, encouraging that the clouds do not much change the overall shape of the 
trajectory in color-color space. This implies that the recovery of surface components is still possible 
given sufficient photometric precision. 

5.2. Modification of the surface reflection spectra by semi-transparent cloud cover 

Our inverse procedure assumes that the observed total albedo spectrum of the Earth can be 
approximated by a linear combination of five distinct spectra corresponding to clouds and the four 
surface components. In reality, the reflection spectrum of a particular surface component covered 
by clouds is not necessarily composed of a sum of the clouds and the surface spectra. 

In order to show the difference more quantitatively, we use rstarGh to compute the reflectivity 
of vegetation entirely covered by clouds of t^u = 10. As plotted in the left panel of Figure [I5l 
the resulting spectrum is best-flt by a sum of clouds, vegetation and soil with fractional areas of 
0.888, 0.279 and 0.135, respectively. The cloud-diluted red-edge feature of the vegetation is partially 
misidentified as soil in this case. The estimated total area exceeds unity so as to compensate for the 
lower reflectivities of vegetation and soil compared to clouds. We repeat the exercise for vegetation 
covered by clouds of = 1; the result is plotted in the right panel of Figure [T5l In this case, the 
spectrum is decomposed as a sum of clouds, vegetation and ocean. 

These misidentifications are the result of our simplified treatment of the template surface 
albedos assumed in the inverse procedure. In reality, the wide range of cloud properties results in a 
variety of cloud spectra (Fig. [8]) that are partially degenerate with the spectra of the other surface 
components. In principle, this issue could be addressed with a more elaborate and realistic model, 
but doing so is not justified without quite high signal-to-noise ratio observations. 
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Fig. 15. — Spectra of vegetation as seen through 100 % cloud cover and the resulting decomposition 
based on model I. Left: the r^d = 10 case. It is best-fit by a combination of % ocean, % snow, 
27.9% vegetation (cross), 13.5% soil (open circle), 88.8% cloud (plus). Right: the Tdd = 1 case. 
It is best-fit by a combination of 10.4% ocean (open circle), % snow, 80.5% vegetation (cross) 
12.0%, % soil, and cloud (plus). 



5.3. Additional time variability 



If an exoplanet is an oblique rotator like the Earth, the overall color of its surface may exhibit 
seasonal variations including those due to autumn leaves and snow cover. Otherwise its surface 
components are taken to be independent of time. Clouds can introduce light-curve fluctuations 
due to their spatial and temporal variations, but the good agreement of the diurnal EPOXI light 
curves and our forward simulations based on a static cloud map (see Section [3]) indicat es that the 



variab il ity due to c l ouds is not important on a 24 hour time scale, as discussed in e.g. iFord et al 



(2001 



Palle et al.l (|2008l ). 



Since we may have to stack the data of exoplanets over a period of a few weeks or more to 
improve th e photon statistics, the cloud vari ations on those longer time scales may well be more 
important (|Ford et al.ll200ll : iPalle et al.ll2008l ). Nevertheless, we may hope that the decomposition 
of the stacked light curves will simply provide us an average measure of each component because 
the light curves are essentially linear combinations of the observing days. 

In turn, the time variability among the light curves folded according t o the diurnal peri od may 
act as a useful probe of cloud pattern and thus weather on an exoplanet. iPalle et al.l (120081 ) found 
that the apparent rotation period is shifted by the net movement of cloud distribution. Also it 
might be possible to use the time variability to separate the contributions from clouds and those 
from surface components by e.g. looking for the day with the lowest reflectivity because it is likely 
to correspond to the least cloudy day; generally surface components have lower reflectivity than 
clouds (Fig. [7]). 
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6. Summary and discussion 

We have presented a methodology to interpret future photometric data of Earth- like exoplanets. 
The methodology is an improvement over the one presented in Paper I in various ways. In particular 
we have incorporated the presence of clouds, utilized the radiative transfer code rstardb for the 
atmosphere and based the reconstruction on a more diverse set of template albedo spectra. In 
addition, we have validated the forward procedure using the EPOXI observations of the Earth's 
light curves as well as using them as an input for the inverse procedure. 

Our main conclusions are summarized below. 



1. Our forward procedure reproduces the EPOXI light curves well using the spatial distribution 
of the cloud cover fraction and optical depth on each observing day. This scheme will be 
useful to predict the light curves of Earth-like exoplanets, potentially detectable with future 
ground-based and space facilities. 

2. While our inverse procedure is approximate and model-dependent, the presence of ocean and 
clouds is routinely detectable from the diurnal light curves of an Earth-twin with S/N> 10. 
However, S/N > 20 is required for a high confidence detection. 

3. For the necessarily limited suite of models considered in Section 14. 5| the systematic uncer- 
tainties associated with the assumptions and input template spectra are similar in amplitude 
to the statistical errors at S/N~10. Thus, systematic uncertainties are not necessarily dom- 
inant for the signal-to-noise levels realistically expected for early exoplanet light-curve data. 
However, it is obviously possible for a terrestrial planet to have important surface components 
without counterparts on the Earth; if so, very serious systematic errors could ensue. 

4. For an Earth-twin exoplanet, S/N > 100 may enable us to detect the presence of components 
other than ocean and clouds in a fairly model-independent fashion. 



As discussed in Section [5l cloud cover adds complexity and uncertainty to the interpreta- 



tion of planetary light curves. Fortun ately, clouds hav e their unic 
reflection spec tra, time variability (e.g 



emis sion (e.g. 



(e.g. iRobinson et al 



Tinetti et al. 
2010l : 



2006a 



Ford et al 



2001 



bl: iKitzmann et al 



de Kok et al 



Palle et al. 



2011 



ue fea tures such as the grey 
20081 ). the effect on thermal 



the anis otropic nature of sc attering 



20111 ). and polarization (e.g. iKaralidi et al.l 120111 ). In fu- 



ture work, it would be desirable to develop techniques for removal of clouds by making use of these 
characteristics in order to extract uncontaminated information on the surface. 

Our identification of surface components depends on the differences in their reflection spectra. 
For example, the presence of ocean is identified as a "dark" components in our inversion proce- 
dure. Although there might be other "dark" components than ocean, this approach will provide a 
useful check for the presence of oceans in combination with other indicators, such as "ocean glint" 
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(IWilliams &: Gaidosll2008l : lOaklev Cashll2009l : iRobinson et al.ll20ld ) or polarimetry (jZugger et al, 



20m . 



In this paper as well as Paper I, we have focused on the case of an Earth-twin, but in order to 
further investigate the potential of photometric light curves for the characterization of terrestrial 
exoplanets, it will be necessary to consider models of planets unlike the Earth in various aspects. 
We plan to investigate such models in future work. 
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A. Pre-processing of the MODIS cloud data 

The MODIS global cloud data available onlinj^ require pre-processing to deal with two of 
their limitations: pixels with invalid cloud cover fraction data and those with invalid cloud optical 
thickness data. 

For the former case, we assign a cloud cover fraction equal to the average value of adjacent 
pixels with valid values, if any. We iterate this procedure to assign values to pixels which initially 



This preprint was prepared with the AAS lATpjX macros v5.2. 
http://ladsweb.nascom.nasa.gov/ 
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have no adjacent pixels with vahd values. Re-binning to 2° x 2° resolution is accomplished by 
averaging. 

In the latter case, pixels may be invalid for one of three reasons: (1) the cloud optical thickness 
is too small to measure; these pixels are assigned a value of 0. (2) the Sun illuminates regions near 
the poles at too low angle of incidence to allow measurements of cloud optical thickness; these 
pixels are assigned values via the same iterative averaging of adjacent pixels used for cloud cover 
fraction (see above). (3) A small fraction of the Earth was not observed by MODIS on the days of 
the EPOXI observations; the resulting gaps were filled by the same iterative procedure. 

Because the radiance is related to the cloud optical thickness in a non-linear way, it is inap- 
propriate to average Tcm in the coarse-bining process. Instead, we assign the value of cloud optical 
thickness in one of the sub-pixels to the whole pixel. 



B. Longitutinal mapping of the surface 



Diurnal variations may be translated into a longit udinal map (ICowan Sz Agolll2008l ) or even a 
2-dimensional map if combined with yearly variations (jKawahara fc Fuiiill20ld ). Here we map the 
diurnal variation of each component reco nstructed from June data (The left panel of Figure [TOj) 
into a 7-slice model (ICowan fc Agolll2008l ) and display the result in Figure [TBI It is apparent that 
some of the major geographical features of the Earth e.g. two oceans, the Sahara desert, and the 
two largest land masses are properly located. 
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VEGETATION SNOW IGBP classification map 




Fig. 16. — The 7-slice longitudinal distribution of each component recovered from the June light 
curves. The indices of the IGBP classification map (the bottom left panel) are — 0: ocean, 1: ev- 
ergreen needleleaf forest, 2: evergreen broadleaf forest, 3: deciduous needleleaf forest, 4: deciduous 
broadleaf forest, 5: mixed forest, 6: closed shrubland, 7: open shrubland, 8: woody savannas, 
9: savannas, 10: grasslands, 11: permanent wetlands, 12: croplands, 13: urban and built-up, 
14: cropland/natual vegetation mosaic, 15: snow and ice, and 16: barren or sparsely vegetated 
( )http://modis-atmos. gsfc. nasa.gov/ECQSYSTEM/index.html ). 



